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Abstract 

We consider a three dimensional, generalized version of the original SPP model for collective 
motion. By extending the factors influencing the ordering, we investigate the case when the move- 
ment of the self-propelled particles (SPP-s) depends on both the velocity and the acceleration of the 
neighboring particles, instead of being determined solely by the former one. By changing the value 
of a weight parameter s determining the relative influence of the velocity and the acceleration terms, 
the system undergoes a kinetic phase transition as a function of a behavioral pattern. Below a critical 
value of s the system exhibits disordered motion, while above it the dynamics resembles that of the 
SPP model. We argue that in nature evolutionary processes can drive the sttategy variable s towards 
the critical point, where information exchange between the units of a system is maximal. 

1 Inroduction 

Collective motion of organisms (e.g. fish schools, bird flocks, bacterial colonies) exhibits a large variety 
of emergent phenomena ||Il|2l[3l|4l|5l|6l|7l[8|. Synchronized motion, symmetrical group formations 
(e.g., V shaped) or swirling patterns emerge in spite of the apparently simple behavioral rules of the 
individual flock members JH [TOl . The self propelled particles (SPP) model was proposed by Vicsek et 
al. inn to describe the onset of ordered motion within a group of self-propelled particles in the presence 
of perturbations. Taking into the effects of fluctuations inevitably present in biological systems was an 
essential generalization of the prior deterministic flocking models such as that of Reynolds |[T2l . The 
original model considers point-like particles moving at constant velocity on a two dimensional surface 
with periodic boundary conditions. The only rule is that, at each time step, every particle approximates, 
with some uncertainty, the average direction of motion of the particles within its neighborhood of radius 
R. This model exhibits spontaneous self-organization; by decreasing the noise parameter, the system 
undergoes a kinetic phase transition from a disordered state to an ordered one where all the particles move 
approximately in the same direction. Due its simplicity and analogy with biological systems comprised 
of many, locally interacting particles, the SPP model soon became a reference model for the flocking 
behavior of organisms ||I3liaiIl[ia[l7l[l8l[l9l|20l|2ll|22l|2l. 

The individual based behavioral rules, determining collective motion, are of particular interest. Im- 
portant elements of behavioral rules are the nature of the perceived information and the affected behav- 
ioral traits. ll24l l25l 26 1. A frequent assumption in models is that the information, perceived by the 
particles, is restricted to the velocity of their neighbors. The interaction range is usually defined by met- 
ric distances, but Ballerini et al. 1241 recently showed that topological distance is the one determining the 
flocking of starlings. The assumption of reflecting on the momentary velocity only may not be enough 
for adequately describing a number of biologically relevant situations. We expect that the behavior of 
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the SPP model is significantly extended if we also incorporate a term corresponding to memory on short 
time scales. This can be achieved by introducing an acceleration term into the equations (rules). This 
is equivalent to separating time scales by assuming that the particles differentiate between two kind of 
information, the first being their actual velocity, the second one corresponding to recent change in their 
direction of motion. For example, in the case of birds, reacting to acceleration may mean that birds can 
give signals to their neigbors about their intended changing of their flight direction by quickly modifying 
their velocity. 

2 Model 

In the three-dimensional, scalar noise (SNM) version ||27| of the original SPP model ifTTll the particles are 
assumed to move with a constant velocity u, and their positions are updated simultaneously according to 

Xi(t + At) =Xi(t) +Vi(t)At, (1) 

where Xj and v, are position and velocity of particle i, respectively. The time increment is set to be 
At = 1. Each particle is assumed to move, with some uncertainty, in the average direction of all the 
particles within a fixed neighborhood of radius R = 1. Hence, the new velocity is given by 

Vi(t + At) = 1/ • M{e, e) • N ((v(t)),,^) , (2) 

where u is the absolute value of velocity, A^(e, ^) is a rotational tensor representing a random perturba- 
tion, (v(t))j /J denotes the average velocity of all particles around particle i within radius R including 
particle i itself, and N(u) = u/|u|. 7W(e,^) performs a rotation of angle ^ around a vector e; ^ is a 
uniform random value in the interval [— r/vr, rjir], whereas e is a random unit vector chosen uniformly 
from the set of vectors perpendicular to N ((v(t))i^R). The order in which ([T|) and Q are calculated has 
some quantitative effects on the results (see later). 

Here we introduce the acceleration coupled self-propelled particles model (AC-SPP) being a modi- 
fied version of SNM, in which the velocity vector Vi(t + At) is a function of both the velocity v(t) and 
the acceleration a(t) = (v(t) — v(t — At))/At of the neighboring particles. Then equation [2] becomes 

vi(t + At) = i^- M{e, • N (s • {v{t))i,R + (1 - s) • (a(t)At)i,/j) , (3) 

where s G (0, 1] is a so-called strategy parameter, expressing the relative influence of the acceleration and 
velocity tags on the velocity vector of the focal particle. Initially, we have N = pL^ randomly distributed 
particles, where L and p stand for box size and particle number density, respectively. The bounding 
box has periodic boundary conditions. The velocity parameter used in the simulations is = 0.1, 
corresponding to the low velocity regime [2T|. 

Taking into account the acceleration in a separate term has various possible motivations. As for the 
contribution of the z-th particle, it corresponds to a memory effect: a given particle, if it has no neighbors, 
has a tendency to keep on turning as it did in the previous time step. Perhaps more importantly, the aver- 
age turning rate of the neighbors has now a separated effect on the turning rate of the i-th particle. In the 
s close to 1 limit the AC-SPP model is very much like the original SPP, while for s ^ 1 the acceleration 
term dominates and instantaneous turning of the neighbors has a strong effect on the trajectory of the i-th 
particle. In this way not only the state (velocity), but a kind of behavioral pattern (turning) can be taken 
into account. We associate, for example, with such a behavioral pattern a short turning period of a bird 
flying in a Hock and giving sign for its neighbors of its intended changing of the direction of flight. 

We characterize the collective motion of particles by the average velocity of all particles 99, defined 
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with denoting the number of particles in the system. This order parameter can take any value in the 
range [0, 1], and expresses the tendency of particles to move in the same direction. If the particles move 
randomly 99 = whereas if every particle moves in the same direction if = 1. (ip) was obtained by 
averaging over 12 individuals runs, each recorded after a relaxation time to stationarity T^eiax and with 
time averaging over an additional Tavr time steps. 

3 Results 

At first we investigated the dynamics of the system at different values of the strategy variable s at fixed 
density and noise values: p = 0.16, 77 = 1/9. The density was choosen to be high enough to get flocking, 
but allowing moderate CPU times per run. The noise value was a typical value resulting in an ordered 
state for the original SPP model {s = 1). We found that by increasing s the system undergoes a phase 
transition; below a critical value, Sc, the ordering, expressed by tp, is negligible, while above Sc the level 
of order increases rapidly as a function of s (Figure [T]). This is a novel type of phase transition since it 
corresponds to a phase change due to a change in the relative strength of a behavioral pattern. 
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Figure 1: Average velocity as a function of the strategy variable s. The data points were obtained 
by averaging the results of 12 simulations, using the parameters L = 100, p = 0.16, v = 1/9, with a 
relaxation time of Treiax = 10000 time steps and averaging over an additional Tavr = 10000 time steps. 

In order to determine the nature of the phase transition we also calculated the probability density 
function (PDF) of ip and the Binder cumulant G at different strategy values. The Binder cumulant, 
defined as G = 1 — {cp^} /3{ip'^)'^, measures the fluctuations of the order parameter and can be used to 
distinguish between first- and second order phase transitions [31 ]. In case of a first order phase transition 
G exhibits a characteristic minimum, whereas in case of a second order transition this sharp minimum 
is absent. In our case, the PDF was unimodal, and G did not have a sharp minimum around Sc, both 
indicating a second order phase transition (Figures [2] and [3]l. G was monotonously increasing, as it was 
observed already in a three dimensional SNM with continuous phase transition [271. 

Consequently, near the critical point, the order parameter obeys the scaling relation 

^r^[s-s,{p,r^)f, (4) 

where (3 is the critical exponent. The critical values Sc and /3 were determined by plotting log{tp) as 
a function of log[{s — Sc)/sc] (Figure |4]). Sc was obtained by finding the value where the plot was the 
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Figure 2: Probability density function (PDF) of the order parameter (/?. Parameters as in the previous 
plot. The one-humped distributions in a range of s values close to the critical one (Sc ~ 0.59) suggest a 
second order phase transition. 
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Figure 3: The Binder cumulant G as a function of the strategy variable s. The dotted lines at 2/3 and 
4/9 indicate the theoretical value in the case of the ordered and totally disordered states, respectively. 
Parameters as in the previous plots. The absence of a sharp minimum indicates a second order phase 
transition. 
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straightest in the relevant region, whereas the critical exponent is equal to the slope of the fitted line. We 
obtained Sc = 0.590 ± 0.002 and /3 = 0.35 ± 0.05. It should be noted, that the value /3 = 0.35 is far from 
0.5 (corresponding to mean- field or bifurcation) indicating the true non-trivial nature of the transition we 
describe. 
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Figure 4: Average velocity values as a function of the distance from the critical point. The fitted line, 
with Sc = 0.590 ± 0.002, had a slope of /? = 0.34 it 0.05. Parameters as in the previous plots. 

The critical value of the transition depends on both the density and the noise parameters. Figure |5] 
shows that Sc is decreasing with both increasing density and noise values. 

The order variable </?, by itself, provides a poor description of the two states, hence we also calculated 
other statistics. The sinuosity of the particle trajectories was expressed by their average curvature, defined 
as A = X^AT ^"'^^p'^ ■ Two other statistics, ip and fi, measure the information exchange between the 
particles, ip is defined as the average number of different particles encountered, i.e., being within a 
distance R, by a focal particle during a given time interval. Although this quantity depends on the 
chosen time interval, its monotonicity as a function of s is independent of it. fi was used to evaluate the 
speed of information propagation, as follows. Initially one percent of the particles held the information. 
The information was transmitted between particles via encounters between information holders and other 
particles. This way sooner or later every particle became information holder, fi was defined as the time 
needed for at least 90 percent of the particles becoming information holders. 

All the three statistics (A), {^p) and (/i) were obtained by averaging over 12 individual runs, each 
with a relaxation time of T^eiax = 10000 and averaging time of Tavr = 10000. (^Z") is an average value 
for 100 randomly chosen particles. 

The curvature of the trajectories decreases with the strategy variable (Figure [6]l. Large curvature at 
low s values indicates in this case, that particles move in small circles (Figure [7]). It is because the large 
influence of the acceleration term results in continuous turning. This turning is synchronized among 
neighbors, i.e., their acceleration and velocity vectors become the same, resulting in a particle cloud 
consisting of separated groups of particles, each containing circling particles. By increasing s the radius 
of these circles increases, until the circling groups overlap and start to interact with each other. At a 
critical point the circles overlap so much, that neither their position, nor their composition remains the 
same; in other words the circling groups lose their identity and the particles start to move sinuously. At 
large s values the movement becomes ordered; all particles tend to move in the same direction, similar 
to the ordered phase of the SNM model at small velocities. 

The dynamics of the system is well reflected in the information propagation (Figure|8]). Both ip and 
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Figure 5: Border surface between the ordered and disordered phases in the s, p, rj parameter space. 
The ordered (((/?) > 0.01) and disordered {{ip) < 0.01) states are indicated by heavy and light dots, 
respectively. Squares show the critical noise values r]c{s, p) obtained from the plot of log{ip) as a function 
of log[{rjc — Tj) /rj^ for each s and p pair. L = 100. 




Figure 6: The curvature of trajectories A as a function of the strategy variable s. The dotted line at 
Sc = 0.590 indicates the position of the critical point for ordering. Parameters as in the previous plots. 
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Figure 7: Positional data and spatial trajectories of particles at different strategy values. Subfigures 
show the typical behavior (a, d) below, (b, e) near and (c, f) above the critical point, respectively, (a-c) 
Positional data of the particles shown as birds, (d-f) Each curve shows the trajectory of a particle over 
60 time steps after reaching steady state (t = 20000). Different shades of gray indicate the time past, 
with darker tones denoting more recent positions, (a, d) Below the critical point, at s = 0.5 the particles 
move in circles, (b, e) Near the critical point, at s = 0.53, the particles move sinuously, (c, f) Above the 
critical point, at s = 0.9, the particles move in the same direction. L = 20, t] = 1/Q, other parameters as 
in the previous plots. 
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Figure 8: Information exchange between particles (ip) as a function of the strategy variable s. The curve 
for /i was very similar (not shown). Parameters as in Figure [T] 
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have low values at small s and have a maximum value around Sc (see also Ref. ISU). This result held for 
all density and noise parameter values investigated. The curves of tp and fj, were very similar, indicating 
that both are proper measures of information propagation. 
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Figure 9: The order parameter (ip, x marks), the Binder cumulant (G, triangles), the information exchange 
between particles (ip, circles) and the average curvature of trajectories (A, squares) as a function of the 
strategy variable s in case of the AC-SPP with SVA. Compared to the OVA in case of the SVA the critical 
strategy value is at lower s, but the behavior of the ip, G and A functions are very similar. The information 
exchange between particles has its maximum value near the critical point as well, however, its variability 
as a function of s is much stronger than for OVA. Parameters as in Figure [T] 

In the original SPP model (TV] (here, following the notation of Huepe and Aldana, Ref {2T\, called 
as the Original Vicsek Algorithm (OVA)) the positions of the particles at t + At are depending on 
two previous time-steps, t and t — At. In the literature some authors have implemented the model in 
a slightly different way |[20l l28l [29l l30l . where the order of the position and the velocity update are 
changed. Huepe and Aldana (22) refer to this as the Standard Vicsek Algorithm (SVA) and they report 
that the local density is different in the OVA and the SVA, while in both cases the average number of 
interacting neighbors is unreasonably high because of the lack of a repulsive effect. By analysing the 
AC-SPP model with SVA-like updating rule (Figure |9]l, we find that the behavior of the order parameter, 
the Binder cumulant and the average curvature are very similar to those in the OVA, but the critical s 
value is lower. The information exchange rate is, however, rather different and the maximum value of tp 
is much higher in case of the SVA. The maximum value of the information exchange compared to the 
value at s = 1 in case of the OVA is V^mL^ / V^fJi^ = 1-4, while in case of the SVA ip^^ /V'f Ji^ = 10.1. 
At s = 1, where the AC-SPP is congruous to the original SPP model, the rate of information exchange 
is similar in the OVA and the SVA. 



4 Conclusions 

In conclusion, we investigated the statistical properties of a three-dimensional self-driven particle system 
(AC-SPP), designed to be an improved model for the collective motion of living beings and possibly 
nonliving units (robots). The ordering of particles exhibited a second-order phase transition as a function 
of the control parameter corresponding to a behavioral strategy in our case. 

We found that the information exchange between particles was maximal at the critical point. Due 
to the important role of information exchange in animal societies, this might indicate that the critical 
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point corresponds to an optimal behavioral strategy. In a more general context, this result implies that 
biological evolution may drive individual traits corresponding to critical values. However, because of the 
individual based optimahty requirement, this possibility needs to be investigated within an evolutionary 
game theoretical framework. 
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